Bloch oscillations of Bose-Einstein condensates: Breakdown and revival 
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We investigate the dynamics of Bose-Einstein condensates (BEC) in a tilted one-dimensional 
periodic lattice within the mean-field (Gross-Pitaevskii) description. Unlike in the linear case 
the Bloch oscillations decay because of nonlinear dephasing. Pronounced revival phenomena are 
observed. These are analyzed in detail in terms of a simple integrable model constructed by an 
expansion in Wannier-Stark resonance states. We also briefly discuss the pulsed output of such 
systems for stronger static fields. 
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I. INTRODUCTION 

Despite its apparent simplicity, the dynamics of quan- 
tum particles in periodic structures is full of surprises, 
even in the one-dimensional case. For almost a cen- 
tury Bloch waves have been known which are delocalized 
states in a lattice leading to transport. If an additional 
static field F is introduced, these states become localized 
and counterintuitively transport is dramatically reduced. 
Instead an oscillatory motion is found, the famous Bloch 
oscillations. These have a (Bloch) frequency u>b — Fd/h 
where d is the lattice constant and they extend over a spa- 
tial interval A/F where A is the width of the first Bloch 
band. During the last decade these Bloch oscillations 
have been experimentally observed which triggered a re- 
newed theoretical interest (for recent reviews see 0,0,|3|). 

For stronger fields decay has to be taken into account. 
So this simple picture must be replaced by introducing 
couplings to higher bands or, alternatively, by a descrip- 
tion in terms of Wannier-Stark resonances. More details 
can be found in the review article 

One of the most interesting systems for exploring the 
dynamics described above are cold atoms in optical lat- 
tices, because here the notorious difficulties met in solid 
state systems (where, in fact, the Bloch oscillations have 
been observed for the first time Q) are absent or, at least, 
can be made very small. Kasevich did one of the first ex- 
periments with a Bose-Einstein condensate (BEC) of 
Rubidium atoms in a optical lattice with gravity acting 
as the static field. They could observe a pulsed coherent 
output of atoms. 

The atoms in a BEC scatter off each other, which of- 
fers the opportunity to study the influence of the atomic 
interaction on the dynamics. In a good approximation, 
the dynamics can be described by the one-dimensional 
Gross-Pitaevskii equation (CPE) |jj 



2M x w 



Fx 



^ (1) 



where M is the atomic mass, g is the interaction strength, 
and V(x) = V(x + d) is the periodic lattice potential. 

This nonlinear system shows basically all the features 
found in the analysis of the linear equation, such as Bloch 
oscillations of the condensate H,@. In addition, the non- 
linearity introduces new effects, such as soliton-like mo- 
tion, nonlinear Zener tunneling [lfj, and "classically" 
chaotic dynamics 0, 0, 0] . This system has been an- 
alyzed with various methods, see e.g. (l5l ITS Il7| . 

In the present article, we will focus on Bloch oscilla- 
tions and analyze the influence of the nonlinearity. The 
paper is organized as follows: In section [H] we present 
results from a numerical solution of the CPE and show 
nonlinear Bloch oscillations for relatively weak fields and 
different strength and signs of the nonlinear interaction. 
Section IIIII introduces our main tool, a discrete repre- 
sentation by an expansion in Wannier-Stark resonance 
states and derives approximate results based on this ap- 
proach. These results are used to analyze the dynamical 
behavior of Bloch oscillations in section IIVI Finally we 
discuss the modification of the coherent pulse output of a 
Bloch oscillating condensate for stronger fields in section 
Ivl The paper closes with some concluding remarks. 



II. NUMERICAL STUDY OF NONLINEAR 
BLOCH OSCILLATIONS 

Due to the nonlinearity of the GPE QJ analytical stud- 
ies are difficult, so numerical simulations are helpful in 
guiding theoretical investigations. 

In all numerical studies we will use a cosine potential 
V(x) — Vo cos(2-Kx/d). We furthermore use scaled units 
with d — 2tt, Vq = M = 1 (see [4| for more details). 
It is worth noting that the scaled interaction strength is 
inversely proportional to the depth of the potential. 

The GPE then reads 



ihdtip = 



h 2 

— dl + cos x + Fx + g \ip\' 



(2) 
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and we will use the value h — 3.3806 for the scaled 
Planck constant adapted to the experiment of the Ka- 
sevich group (see also 0])- The nonlinearity param- 
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eter is of the order g < 1 in this experiment. Here we 
will extend the analysis, however, to much stronger non- 
linearities up to \g\ = 10. This regime could be reached 
experimentally by increasing the transverse confinement 
or decreasing the depth of the optical lattice. 

Here we are mainly interested in the dynamics of Bloch 
oscillations and therefore use a weak field (F = 0.005) 
and initial states populating almost exclusively the low- 
est Bloch band. In this case the decay is negligible. In 
the linear case the band gap between the lowest and the 
next higher Bloch band for the field free case is 5 = 0.998 
and the probability for Zener tunneling is ~ 10~ 12 . This 
does not change noticeably in the nonlinear case as has 
been checked numerically, showing that nonlinear Zener 
tunneling [ToL Hll | does not play a role for the given pa- 
rameters. Furthermore note that the dynamical effects 
which are observed in the following examples (e.g. the 
damping of Bloch oscillations) are essentially captured 
by a simple integrable model which is introduced in sec- 
tion [ffi] This should be compared to the breakdown of 
Bloch oscillations due to dynamical instability discussed 
in [HEHl. 



Bloch oscillation for the linear case 3 = 0. As an initial 
state, we use a Gaussian wave packet 




FIG. 1: The squared modulus |^(x,f)| 2 of the wave function 
for g = shows the familiar Bloch oscillations. 
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FIG. 2: Expectation values of the position (x) t and width 
Axt for the wave function shown in figure 

Let us start our discussion with a brief look at the 



ip(x,t= 0) = 



1 



(2 7 r) 1 /4aV2 



-(x-x ) 2 /4a 2 



(3) 



with width a = AOtt which is then projected onto the 
lowest Bloch band and normalized to unity. For the time 
propagation a split-operator method .21] is used which 
can also be applied to the nonlinear case. In figure Q we 
observe the familiar Bloch oscillation with a large am- 
plitude because of the weak field. Let us recall that the 
region over which the Bloch oscillation extends can be 
estimated as A/F sa 200 « 32 • 2ir within the tilted band 
picture, where A = 0.9994 is the width of the dispersion 
relation E(k) in the field free case. The numerical re- 
sults confirm this estimate as the top of figure [21 shows. 
As expected for such an initially wide distribution in co- 
ordinate space, the width of the wave packet remains 
practically constant varying periodically with a relative 
amplitude of about 10~ 3 (figure |21 bottom) . 




FIG. 3: Same as figure^ however for a nonlinear interaction 
g = +5 (top) and g = — 5 (bottom). 

Let us now discuss the influence of a nonlinearity fix- 
ing for the moment the nonlinear parameter at g = 5 
(repulsive interaction) and g — —5 (attractive interac- 
tion). From figure |3 one can observe that the the Bloch 
oscillations continue to exist, at least for the short times 
up to t w 10 Tg shown in the figure. In addition to 
the well known localization of \i/>(x, t)\ 2 in the minima of 



3 



the cosine-potential, one observes a further filamentation 
which is particularly pronounced for an attractive inter- 
action. As shown in figure 0] the amplitude of the oscil- 
lation decreases and the oscillation in the width strongly 
increases. 
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FIG. 4: Expectation values of the position (x) t and width 
Axt for the wave function shown in figure [3] for a repulsive 
nonlinearity g = +5. 
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FIG. 5: Same as figure |1J however for for an attractive non- 
linearity g = — 5. 

Also shown in figures 0] and is the time dependence 
of the width Ax t of the wave packet. In sharp contrast 
to the tiny oscillations of the width in the linear case (see 
figure |5J| we find here very pronounced oscillations which 
are rapidly growing (as already described by Holthaus 
|l5|h Such a phenomenon is known as breathing and 
is exhibited in the linear system by wave packets which 
are initially strongly localized in coordinate space pll2^. 
Note that the oscillations of the width Axt for a repulsive 
and attractive nonlinearity are opposite to each other. 

For a stronger nonlinearity g — 10, as illustrated in 
figures and the Bloch oscillations are damped more 
strongly. However, the oscillation docs not fade com- 
pletely but shows a revival with a smaller amplitude af- 
ter a shrinking to approximately two lattice periods. A 
corresponding behavior is observed for the width, where 
the breathing amplitude of the wave function first grows 



fast up to a time of about eight Bloch periods. After 
this time, the width remains limited and oscillates in the 
interval from 31 to 35 lattice periods. 
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FIG. 6: Expectation value (x) t of the position and width 
Axt for a repulsive nonlinearity g = +10. The inset shows a 
magnification of the time interval between 8 and 16 Tb ■ 




FIG. 7: Width Ax t of the wave packet shown in figure|S|for 
a repulsive nonlinearity g = +10. 

Furthermore, the oscillation of (x)t shows phase jumps 
that can be seen, e.g. in the inset of figure|f)]at t 14 Tb. 
A similar behavior has also been described in [l5T |. This 
phase jump coincides with a minimum in the amplitude. 
These phenomena can be understood in terms of an ex- 
pansion in Wannier-Stark basis functions as explained in 
the next section. 



III. WANNIER-STARK BASIS SET 
EXPANSION 

An alternative approach to a direct numerical integra- 
tion of the GPE is an expansion in an adequate discrete 
basis as for example the ground states of single potential 
wells In this work we adopt a different approach, 

following |l3tll^|. and expand the wave function into the 
resonance eigenstates of the linear system, the so-called 
Wannier-Stark states 9 a ,n{ x ) which are eigenstates of 
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the linear Hamiltonian Hq: 



(4) 



where a is the ladder index and n is the site index. The 
energies form the Wannier-Stark ladder 



£a.n = £a,Q + 2'KFn. 



(5) 



The Wannier-Stark states extend over several periods of 
the potential (see remark p3 | and review £| for more in- 
formation). This approach has proven to be extremely 
convenient to describe the dynamics in tilted optical 
lattices in the linear case, especially for higher field 
strengths jlll^. 

Up to section E] we will restrict the discussion to 
small field strengths F. Then one can neglect decay and 
Landau-Zener tunneling and use the lowest ladder a = 
only, henceforth the index a is omitted. Also neglecting 
decay, the imaginary part of the energy £ n is set to zero. 
Plugging the expansion ip(x,t) = ^ m c m {t)^ m {x) into 
the GPE QJ leads to a set of coupled ordinary differen- 
tial equations: 

ihy^Cm^m = y^(£p + 2-KFm) Crn^m 
m m 

+s5^cjqc ro *j*,* t „. (6) 

klm 

The energy £q only leads to a global phase factor and 
hence is omitted in the following. The Wannier-Stark 
states \P n are orthogonal to their left eigenstates iff^n for 
m 7^ n. Nevertheless, since we neglect the resonance 
properties of the system we can identify left and right 
eigenvectors, i.e. assume that H a is hermitian. So mul- 
tiplying equation 10 by and integrating yields 



ihc n — 2kF nc n + g 



Xklm c k c l C mi (7) 



klm 



with the coupling tensor 



Xklm = J *;(»)*k(it)*l(l)*m(l)dl, 



(8) 



which is symmetric under the exchange of its first and last 
two indices. Due to the discrete translational invariance 
of the Wannier-Stark states ^> n (x) — ^>q(x — 2im) one 
finds 



ihc n = 2-KFnCn + .9 ^ Xki m c* k+n ci +n c m+n , 



(9) 



kin 



— Xklr, 



defining Xklr, 

Though not suited for direct numerical calculations be- 
cause of the triple infinite sum, equation |JSJ provides a 
basis for further approximations. In the following we will 
reduce it to a simple integrable model, which neverthe- 
less captures important features of the dynamics. To this 



end we decompose the coefficients c„ into phase and am- 
plitude 



l<Pn 



(10) 



The imaginary parts of the coupling tensor Xkim are neg- 
ligible and so one arrives at the coupled equations 



tup n = —2'KFn 



9Pn 2_^Xklr. 

kirn 



Pk+nPl+nPm+n 
rn 



h() n = 2gp 2 n ^2xklr 



x cos((pi+ - fk+n - <Pn), 

1/2 



1/2 



(11) 



kin I 



Pk+nPl+nPm+n 
Pn 



X Sm(fi +n + Lpm +n - Iflk+n - iPn)- (12) 

If the initial state is broad, populating about 20 wells, 
the amplitudes p n (t = 0) are small. Because of p n ^ p 2 nl 
this implies that the amplitudes p n change only slowly in 
time compared to the phases ip n and can be assumed to 
be constant. 

Furthermore we reduce the expression for (p n to the 
most important contributions. Numerically examining 
the Xkim shows that the dominating terms are Xooo, 
XkkQ = Xkok and Xokk, which is not unexpected consid- 
ering equation ©. It can also be argued (and verified 
numerically) that the terms in equation (|11|) which have 
a nonzero argument of the cosine have little importance, 
as their contributions average out. This leaves the terms 
including xooo and Xkko = Xkok and finally leads to: 



fapn = —2'KFn - g^nPn 



with 



In 



Xooo 



fc#0 



XkOk- 



Pn+k 



These equations (|13f) are integrated to 

Pn(t) = p n (0) , ip n (t) = U) n t 

with 



-2-kFu 



91nPn- 



(13) 



(14) 



(15) 



(16) 



Note that this solution is exact for g — 0. Numerical 
calculations show that one can safely neglect the depen- 
dence of 7„ on the index n and set 7„ ~ 7. For the 
given parameters equation i fTljl yields 7„ > 70 = 0.278. 
However, the best fit with the results from a wave packet 
propagation are obtained for 7 = 0.15. 

This admittedly quite crude approximation shows very 
good agreement with an exact numerical solution. In 
figures |H1 and [§] the approximation l|15|) is compared with 
the results obtained by a wave packet propagation using 
the split-operator method [2]]]. A normalized Gaussian 
initial distribution with coefficients 



-n 2 /4a 2 



a = a /2k = 20, 



(17) 
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is used which closely resembles a Gaussian wave packet 
projected onto the lowest Bloch band in configuration 
space. The dynamics for a moderate nonlinearity .<? = 1 
is well described by equation H15|) , only the growth of the 
width is somewhat underestimated. For g — 10 the ap- 
proximation l|15|) becomes less accurate. In particular it 
overestimates the revival of the Bloch oscillation and un- 
derestimates the growth of the width of the wave packet. 
However it still captures the important features at least 
qualitatively: the decay and revival of the oscillations, 
the phase jump around t = 14 Tg and the breathing of 
the wave function. 




FIG. 8: Bloch Oscillations: Expectation value of position 
{x)t of the wave packet for a moderate nonlinearity g = 1. 
The propagation was done with the split-operator method 
(dashed line) resp. with approximation 1151 (solid line). 
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FIG. 9: As figure [S] but for a strong nonlinearity g = 10. 

The systematic growth of the width of the wave packet 
is mainly due to a broadening of the amplitude distribu- 
tion p n and therefore clearly not included in approxima- 
tion l|15|) . To discuss this effect we briefly reintroduce the 
time-dependence of the p n . Again we reduce the triple 



sum to keep the calculations feasible. Note that the terms 
are oscillating due to the sine. Using approximation 115|) 
we see that for k = I + m the terms proportional to F in 
the argument of the sine cancel and hence the sine oscil- 
lates slowest. Thus we approximate the dynamics of the 
amplitudes p n by 

/ \ V 2 

hp n 2gp n Xl+m,l,m I 5 I 

l,m \ Pn J 

x s\n{tp l+n + ipm+n - <Pl+m+n ~ <Pn), (18) 

where the sum can be truncated at \l\, \m\ = 30. Equa- 
tion (|18|l for p and equation (|13J) for (p are now solved 
numerically with 7 = 0.15 and the initial condition (|17|l . 
The results displayed in figure HUI show that this model 
captures the growth of the width of the wave packet. We 
will, however, not go into details here and return to the 
approximation (|15|) to discuss the dynamics of Bloch os- 
cillations. 
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FIG. 10: Bloch Oscillations for a strong nonlinearity g — 
10. The propagation was done with split-operator method 
(dashed line) resp. with approximation 1181 (solid line). 



IV. ANALYSIS OF THE DYNAMICAL 
BEHAVIOR 

Further insight can be provided by a closer look at the 
dynamics of the wave function in momentum space. This 
can be achieved with the approximate time evolution of 
the expansion coefficients c n derived in the previous sec- 
tion. 

First we briefly reconsider the linear case. For g = 
the equations (|15fl reduce to: 

p n (t) = p„(t = 0) and ip n (i) =-2%Fnt/h. (19) 

The Wannier-Stark functions [24| are related by a 
spatial translation ^ n (x) — ^!o(x — 27m). In momentum 
space this reads 

* n (A) = e- i2 ™ fc * (fc) (20) 
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and the time evolution of the wave function in momentum 
space is 

ip(k,t) = y {k)^2yfp^eKp(-i2im(k + Ft/h)) 

~ V (k)C{k + Ft/H), (21) 

neglecting a global phase. Thus the wave function is the 
product of a time-independent function ^o(k) and the 
discrete Fourier transformation C(k) of the amplitudes 
evaluated at the point k + Ft/H. 

The function C{k) is periodic in momentum space: 
C(k + n) = C(k) for n £ Z. Thus the function 
C(k + Ft/H) is periodic in time with the Bloch period 
Tb = H/F. For a broad Gaussian distribution of the am- 
plitudes p n the discrete Fourier transform C(k) is a comb 
function with narrow peaks at k = n. 

So one arrives at a simple view of the dynamics in mo- 
mentum space: The comb function C(fc) moves uniformly 
under the envelope ^o(k), as illustrated in figure ITT1 In 
coordinate space this periodic motion appears as a Bloch 
oscillation ID. 
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FIG. 11: Bloch oscillations in momentum space. The wave 
function \ip(k, t)\ shown for t = (solid line) moves uniformly 
under the envelope of the Wannier-Stark function [\Po(fc)| 
(dashed line). 

In the nonlinear case one has to evaluate 
t/>(M) = *o(fc) ^2 VPn ex P (~i ( 27r nk - <p n {t))) 

n 

= * (k)C(k,t). (22) 

instead of equation (|21|l . The time evolution of the phases 
is approximated according to equation (|15|) : 

<Pn{t) = Lo n t with fujj n — —2irFn — gjp n , (23) 

and the amplitudes are assumed to be Gaussian 



Pn{t) 



-n 2 /2a 2 



(24) 



As in the linear case, the static field term —2^Fn in equa- 
tion (|23[l for the frequency leads to a uniform motion of 
the function C{k, t) in momentum space. The nonlinear 



term — g^p n leads to a dephasing of the coefficients c n and 
broadens the Fourier transform C(k,t). This dephasing 
causes a damping of the Bloch oscillations in coordinate 
space. 
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FIG. 12: Time evolution of the function C(k, t) in equation 
for g = 10. The function C is scaled as (5(0, 0) = 1. 

The oscillations of the width Axt, the breathing, can 
also be understood with this approach. In the linear case 
such breathing occurs for wave functions that are initially 
strongly localized in coordinate space and thus have a 
broad momentum distribution. As explained, the nonlin- 
ear term leads to a broadening of the function C(k, t) and 
hence to a broadening of the wave function in momentum 
space and breathing in coordinate space. For even longer 
times the coefficients dephase completely and the oscilla- 
tions in both the position and the width are damped. 

In order to understand the revivals and the phase 
jumps in the oscillations of (x)t, we need to look at the 
time evolution of the function C(k,t). In figure IT21 the 
function \C(k, t)\ for g — 10 is plotted at times t = 0, 
t = 6Tb, t = 12 T B and t = 18T B . The dynamics of 
the expansion coefficients was calculated with equation 
(|15fl resp. (|23|l . One observes that the initially narrow 
peaks are broadened and an oscillatory structure devel- 
ops with two maxima at the edges of the band populated 
by |C(fc,£)|. These maxima eventually merge, leading to 
a revival of the Bloch oscillations. The new maximum af- 
ter the merger is displaced by Afc = 0.5 in comparison to 
the linear case and hence the phase of the Bloch oscilla- 
tions is reversed. This maximum broadens again, leading 
to a periodic breakdown and revival. The phase of the 
Bloch oscillations is reversed after each breakdown and 
the amplitudes of the revivals decrease. However, these 
further revivals are observed only within the approxima- 
tion (|15f) and not in a wave packet propagation. 

Now we consider the time dependence of the expecta- 
tion values of position and width. These quantities can 
be evaluated analytically in the linear case g = using 
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a tight -binding approximation 0, • In this approx- 
imation the expectation value of the position oscillates 
harmonically with the Bloch frequency ojb, 



= x + Acos(ujBt) 



and amplitude 



A = 



-2ir 2 Ap 2 /h 2 



2F 



(25) 



(26) 



where Ap is the width in momentum space and A is the 
bandwidth of the dispersion relation E(k) in the field- 
free case. For a small nonlinearity the broadening of the 
wave function in momentum space due to the nonlinearity 
happens slowly compared to the Bloch oscillations. Thus 
we can assume that (x) t still executes damped harmonic 
oscillations with the amplitude (|26|l . where the damping 
is determined by the slowly increasing momentum width 
Ap t . 

According to equation IL'l'l) we can estimate the mo- 
mentum width Ap t by the width of the peaks of the 
function C(k,t). For a broad distribution of the coef- 
ficients c„, as assumed throughout this paper, the sum 
in equation 12211 can be replaced by an integral: 



C(k,t) 



/p^exp (— i(27rnfc — ip n (t))) dn. (27) 



This expression is valid for \k\ < 0.5, otherwise C(k,t) 
is determined by its periodicity. The amplitudes p n and 
phases ip n (t) are approximated by equation 115(1 . where 
the amplitudes p n are normalized as 



1 



V2~7 



:e -n J /2o 



(28) 



We note that C(k,t) depends on the momentum k only 
through the expression k = k + Ft/h, reflecting the uni- 
form motion in momentum space due to the static field: 



C(k) = 



</2 7 



/•+00 / 2 

' & -n 2 /ia 2 e -i(27rfcn+0o-" ' 2c 



'dn 
( 29 ) 

with ft = -ygt/(V2Tr<rh). 

The integral l|27|) can be evaluated using the stationary 
phase approximation. However, there exists only a finite 
fc-interval for which stationary points exist. For |fc| > \k c \ 
with 



2nae 1 / 2 



(30) 



the integral vanishes in the simple stationary phase ap- 
proximation [27j. For \k\ < \k c \ one obtains 



C(k,t) 



2 \ V4 



V /3 2 



e iC+ 



e iC- 



(31) 



where z± are the two solutions of the equation 

2 



ze 



(32) 



! 1/2 

and the abbreviations C± = —2irak (z± i/z + z± 1 ) were 
used. As an example the function \C(k,t)\ is plotted in 
figureEUfor gt = 90 Tg and 7 = 0.15. 




FIG. 13: Function \C(k,t)\ (23 for gt = 90T B and 7 = 0.15. 
The integral was evaluated numerically (solid line) resp. using 
the stationary phase method (dashed line). 

Estimating the momentum width as Ap w h\k c \ one 
arrives at 



Apt 



hg\ 



(27r) 3 / 2 (7 2 e 1 /2 



t. 



(33) 



Thus the damped Bloch oscillations in coordinate space 
are described by 



A 

{x) t ~x + —exp 



i 2 g 2 t 2 



cos(w B t) (34) 



according to l(25|l . The amplitude decreases exponential!; 



laity 

in 



with —g t in agreement with the estimate given in 

Furthermore we can calculate approximately the time 
up to the first rephasing of the coefficients and thus the 
first revival of the Bloch oscillations. This revival occurs 
if the outer peaks of C(k,t) meet at k = n + 1/2, n € Z, 
as illustrated in figure IT21 Therefore the first rephasing 
and revival occurs for k c = 0.5 which yields a revival time 



_ (2^) 3 /V/ 2 ha 2 
2|7Sl 

For g — 10 and 7 = 0.15 one obtains 



tr 



17 T B 



(35) 



(36) 



in reasonable agreement with the revival of Bloch oscil- 
lations observed numerically for the wave packet propa- 
gation shown in figure 

For very long times the coefficients c„ dephase com- 
pletely. We can therefore estimate the position expec- 
tation value by approximating the wave function as an 
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incoherent sum of the basis states. Assuming that the 
amplitudes p n are constant in time according to equa- 
tion (fH5)l one has 

(x)=J2Pn(*n\x\y n ). (37) 

n 

Using the translational properties of the Wannier-Stark 
states (cf. |4|) one arrives at 

{x) 00 a(9 \x\9o) + 2ir^2np n . (38) 

n 

The amplitudes of the initial state (|17fl are symmet- 
ric around n — and hence this approximation yields 
(xjao w (\E'o|a;|^ / o) = —10.5 • 2ir. This estimation fairly 
agrees with the numerical results displayed in figure|SJ As 
argued above, the systematic growth of the width of the 
wave packet is mainly due to a broadening of the ampli- 
tude distribution p n and hence cannot be explained using 
the simple model discussed here. 



V. STRONG STATIC FIELD AND DECAY 

An expansion into Wannier-Stark resonances is also 
very helpful in order to understand the dynamics and 
decay in strong static fields. In the following we will dis- 
cuss the dynamics for the parameters K — 3.3806 and 
F = 0.0661, corresponding to the experiment of the Ka- 
sevich group jfj . A detailed discussion of this experiment 
in terms of Wannier-Stark resonances, however neglect- 
ing the nonlinearity, can be found in hj. Thus we will 
only briefly discuss the influence of the nonlinearity on 
the pulse shape. 

For a field strength of F = 0.0661 decay cannot be 
neglected any longer. One has to take into account that 
the resonance states eventually diverge exponentially for 
x or k — > —oo. Hence, a wave function of the form ill' 21) 
is not normalizable. Nevertheless, the restriction to the 
ground Wannier-Stark ladder is still sufficient. 

As described in [l^ (see also remark [24j) one can solve 
the problem of normalization by introducing truncated 
resonance states defined by 

*K(k) = Q{k + K)9 n (k). (39) 

The Heaviside-function 0(fc + K) truncates the reso- 
nances at —K. Provided that \K\ is large enough, the 
time evolution of these states is given by 

¥*(fc, t) = 0(A + K + Ft/h)V n (k, t). (40) 

If the support of the initial wave function is bounded in 
momentum space by |fc| < \K\, we can expand it into a 
basis of truncated resonances. The dynamics of this state 
is then given by 

ii{k, t) = e{k + K + Ft/K) tf (fe, t) C(k, t), (41) 



with ^o(k,t) — exp(—i£ot/h)'i'o(k) instead of equation 

For a coherent initial distribution of a sufficient jvidth 
c n ~ exp(— n 2 /(2a) 2 ) with u > 1, the function C(k,t) 
is a comb function in the linear case, leading to a pulsed 
output. The pulse shape given by the function C broad- 
ens and deforms under the influence of the nonlinearity 
as described in the previous section (cf. figure IT2*)l . This 
deformation is directly observable in the pulsed output. 




k 

FIG. 14: Pulsed output for different nonlinearities g = 
0, —5, 5 and 10 (from top to bottom). The wave function 
|V>(M)| 2 is displayed for t = 9.1Tb- 

This is demonstrated in figure ITU for a coherent initial 
distribution c„ ~ exp(— n 2 /4cr 2 ) with a width of a = 
15/2. The time evolution was again calculated using a 
split-operator method. The wave function \ip(k, t)\ 2 is 
plotted at a time t = 9.1Tb for four different values of 
the nonlinearity g — 0, — 5, 5 and 10. The broadening 
of the peaks with increasing \g\ and the characteristic 
double-peak structure can clearly be seen. 

The resulting wave function in coordinate space is a 
sequence of pulses at the points 

x = x -^(t + jT B ) 2 , (42) 

with xq = Eg /F and j E Z. These pulses are accelerated 
just like classical particles in a static field, as observed 
in the experiment 0. The pulse shape is approximately 
described by the discrete Fourier transform of C fl8j| . 
Thus one also finds a characteristic deformation of the 
pulses in coordinate space. 

VI. CONCLUSIONS 

In this article we first investigated Bloch oscillations 
of BECs by numerical solutions of the Gross-Pitaevskii 
equation (GPE) and demonstrated a revival of Bloch 
oscillations after an initial breakdown. These findings 
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have been further analyzed via discretising the GPE in 
a Wannier-Stark basis set expansion. Using these reso- 
nance states one can easily compare the linear and nonlin- 
ear case. This comparison leads to a better understand- 
ing of the nonlinear features of BECs in optical lattices. 
It allows us to derive a simple integrable model Ijl5(l which 
can explain the nonlinear phenomena of breakdown and 
revival of the Bloch oscillations. This approach, unlike 
the tight-binding approximation, works as well for strong 
Stark fields. 

Many interesting questions are left open and deserve 
future studies, as for example the following: (a) The ef- 
fects induced by the nonlinearity for Bloch oscillations 
in two-dimensional lattices, where recently novel effects 
concerning the extreme sensitivity on the field direction 
with respect to the lattice have been found [3, |25| . (b) 
For BECs in tilted optical lattices a classical chaotic be- 



havior has been reported [13| . These interesting findings 
deserve further studies, for example the correspondence 
of the emergence of chaos with the loss of stability of 
the GPE solutions, (c) It is an entirely open question 
how the nonlinearity influences the behavior of a driven 
Wannier-Stark system, e.g. the stabilizing phenomena 
found for an additional harmonic driving [4], l2f)| . 
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